Gain Drifts and White noise

When turning on the noise we have to make sure that the white noise component of gain fluctuations $\sigma _g $ is lower than the white noise level $\sigma_n$ which is set by the detector NET, $\sim 40 \mu K \sqrt{s} $ (gain drifts would add unphysical noise otherwise ).

Moreover, we added to the overall input signal the CMB monopole being the major contributor to the thermal load on detectors.

We therefore simulate observations scanning just the monopole and injecting gain drifts with $f_{knee} = 30 $ mHz. We then stored in memory the timestream relative to 1 precession period with and without the noise coaddition (see the plot below).

In [28]:
import pylab as pl 
from  numpy.random import normal 
import numpy as np
import healpy as hp 
from scipy.signal import periodogram 
from pylab import cm 
cmap = cm.get_cmap('RdYlBu')
cmap.set_under('w')
cl = pl.loadtxt('/Users/peppe/work/satellite_sims/Class_planck2015_r0_cl_lensed.dat')

Nl =lambda w: (pl.pi /180. /60. * w)**2
wp=7.6 # uK arcmin 
Nl0= Nl(wp)
In [3]:
from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[3]:
In [99]:
fs =40. 
fk =0.03 
NET= 40e-6
wnoise= pl.load('/Users/peppe/work/satellite_sims/drifting_gain/tod_after_wnoise.npy')
wonoise= pl.load('/Users/peppe/work/satellite_sims/drifting_gain/tod_after.npy')
noise = pl.load('/Users/peppe/work/satellite_sims/drifting_gain/tod_before_m_wnoise.npy')
function=lambda  f,wn, fk :wn *(1+  fk/f )
f,wonoise_spec =periodogram(wonoise , fs=fs, scaling='density'  ) ;

f,wnoise_spec =periodogram(wnoise, fs=fs , scaling='density'  )

f, noise_spec =periodogram( noise , fs=fs, scaling='density'  ) ;

t=pl.arange(len(wnoise)) /fs/3600.
pl.figure(figsize=(10,5))
pl.subplot(121) 
#pl.plot(t, wonoise, alpha=.5, label='CMB Monop. w/drifts')
pl.plot(t, wnoise, alpha=.5, label='CMB Monop. + gain drifts + noise ')
pl.plot(t, noise, alpha=.5, label='CMB Monop.  + noise ')

pl.xlabel('t [ h ]', fontsize=15)
pl.ylabel('TOD [ K ]  ', fontsize=15)
pl.legend() 
pl.subplot(122);
pl.ylim([1e-7, 1e-1]); 
pl.xlim([1e-4, 1e1]); 

#pl.title('') ; pl.loglog(f, pl.sqrt(wonoise_spec), alpha=.5 , label='CMB Monop.  +drifts');

pl.loglog(f, pl.sqrt(wnoise_spec), alpha=.5 , label='CMB Monop. + gain drifts +noise');
pl.loglog(f, pl.sqrt(noise_spec), alpha=.5 , label=' CMB Monop. + noise');

pl.loglog(f[1:], function(f[1:] , NET ,fk) , lw=3, alpha=.7,label=r'$f_{k}=30\, mHz , \, NET= 40 \mu K \sqrt{s}$' );
pl.legend() 

pl.grid(True)
pl.ylabel('PSD [ '+r'$K \sqrt{s}$' + ' ]', fontsize=15)
pl.xlabel('f [Hz]', fontsize=15)
Out[99]:
<matplotlib.text.Text at 0x95a784850>
We 've found that $\sigma_g =10^{-3}$ is enough to simulate gain drifts whose fluctuations are the same order as the noise .

Runs w/ $f_k = 30 $ mHz

Inputs: Dipole , CMB T (No noise)

we simulated signal-only simulations including CMB T and Dipole. We didn't include noise in the simulations since the residual are about 3 orders of magnitude lower than the noise level coming from 2 detector pairs. We should rescale the noise level to assess the level of nominal noise wrt the residuals coming from gain drifts.

In [115]:
umap1= hp.read_map(field=range(3) , 
            filename='/Users/peppe/work/satellite_sims/drifting_gain/dipole/uncorr/madam_map_fk0p03.fits', 
                       verbose=False)

lmax=3*hp.get_nside(umap1)
l=pl.arange(2, lmax+1)

fig=pl.figure(figsize=(15,8))
fig.suptitle('gain drifts injected ')

hp.mollview(umap1 [1], unit='K ', title=r'$  Q$', sub=221 ,max=2e-8, min=-2e-8, fig=fig,cmap=cmap  ) 

hp.mollview(umap1[2], unit='K ', title=r'$U$', sub=222 ,max=2e-8, min=-2e-8, fig=fig ,cmap=cmap   )

cmap1= hp.read_map(field=range(3) , 
            filename='/Users/peppe/work/satellite_sims/drifting_gain/dipole/corr/madam_map_fk0p03.fits', 
                       verbose=False)

fig=pl.figure(figsize=(15,8))
fig.suptitle('w/ mitigation')
hp.mollview(cmap1 [1], unit='K ', title=r'$  Q$', sub=221 ,max=2e-8, min=-2e-8, fig=fig,cmap=cmap  ) 

hp.mollview(cmap1[2], unit='K ', title=r'$U$', sub=222 ,max=2e-8, min=-2e-8, fig=fig ,cmap=cmap   )

Fig1: (top) Q and U maps with gain drift injected with a $f_k =30 $ mHz and $\sigma_g =0.1$ %. (bottom) Q and U maps with mitigation to gain drifts applied. Notice that the mitigation is able to reduce the large scale features observed in Q and U maps in the top panel.

Runs w/ $f_k = 150 $ mHz

In [116]:
umap2= hp.read_map(field=range(3) , 
            filename='/Users/peppe/work/satellite_sims/drifting_gain/dipole/uncorr/madam_map_fk0p15.fits', 
                       verbose=False)

fig=pl.figure(figsize=(15,8))
fig.suptitle('gain drifts injected ')
hp.mollview(umap2 [1], unit='K ', title=r'$  Q$', sub=221 ,max=2e-8, min=-2e-8, fig=fig,cmap=cmap  ) 

hp.mollview(umap2[2], unit='K ', title=r'$U$', sub=222 ,max=2e-8, min=-2e-8, fig=fig ,cmap=cmap   )

cmap2= hp.read_map(field=range(3) , 
            filename='/Users/peppe/work/satellite_sims/drifting_gain/dipole/corr/madam_map_fk0p15.fits', 
                       verbose=False)

fig=pl.figure(figsize=(15,8))
fig.suptitle('w/ mitigation')
hp.mollview(cmap2 [1], unit='K ', title=r'$  Q$', sub=221 ,max=2e-8, min=-2e-8, fig=fig,cmap=cmap  ) 

hp.mollview(cmap2[2], unit='K ', title=r'$U$', sub=222 ,max=2e-8, min=-2e-8, fig=fig ,cmap=cmap   )

Fig2: (top) Q and U maps from gain drifting with a $f_k=150$ mHz and $\sigma_g= 0.1$ % . (bottom) Q and U maps with the mitigation applied. The leakage amplitude is larger as compared to the case with $30$ mHz. However, post -mitigation residual are comparable as in Fig1.

BB- Power Spectra

In [70]:
cl1=hp.anafast(umap1 , lmax=lmax)
clc1=hp.anafast(cmap1 , lmax=lmax)
cl2=hp.anafast(umap2 , lmax=lmax)
clc2=hp.anafast(cmap2 , lmax=lmax)
In [82]:
import scipy as sp
from scipy import interpolate
from scipy import optimize
from scipy import integrate
from iminuit import  Minuit
import numpy as np 
def get_delta_r(inputcl, lmax=200, wp = 7.6, fsky=1.,b=1.e-2, n_steps=1024):
        lmin=2
        
        ell=pl.arange(lmin, lmax)
        tenscl =interpolate.interp1d (ell, 
                 pl.loadtxt('/Users/peppe/work/satellite_sims/Class_planck2015_r1_cl_unlensed.dat')[:,3][lmin:lmax])
        lenscl =interpolate.interp1d (ell, 
                pl.loadtxt('/Users/peppe/work/satellite_sims/Class_planck2015_r0_cl_lensed.dat')[:,3][lmin:lmax ])
        nl_w=interpolate.interp1d (ell, inputcl[lmin:lmax]   )
        clobs   =lambda  r, l : lenscl(l) +  nl_w(l)
        cltheor =lambda r,l:  r *tenscl(l)+    lenscl(l)   
        logP    =lambda r, r_inp,l: -   0.5*(2.*l +1 )*(clobs(r_inp,l)/cltheor(r,l)  +pl.log(cltheor(r,l)) -(2. *l -1)/(2.*l+1) *pl.log(clobs(r_inp,l)))
        func=lambda x: - pl.sum( logP(x, 0, ell)  )
         
        #m=Minuit(func, x=1e-3 , error_x=0.0001, limit_x=(0.,1e-2), errordef=1e-8 , print_level=0)
        m=Minuit(func, x=1e-3 , error_x=0.0001, limit_x=(0.,b), errordef=1e-8 , print_level=0)
        
        m.migrad()
        rmax= m.values['x']
        L       =lambda r: pl.exp ((logP(r,0,ell)   ).sum() - ( logP(rmax , 0, ell)  ).sum() )
        Lmap= map(L, pl.linspace(0,b, n_steps ) )/L(rmax)
        a=0
        
        integr_steps=n_steps
        Normalization= integrate.quad(L, a=a, b=b , limit=1000, epsrel=1e-8)[0]
        R=0
        dr= (b-a)/integr_steps

        deltar=a
        it_counter=0
        while R<0.6801 and deltar<b :
                deltar+=dr
                R=  integrate.quad(L, a=a , b=deltar, limit=1000, epsrel=1e-8)[0] / Normalization
                it_counter+=1

        deltar/=pl.sqrt(fsky)
        #print "Delta r= %.2f x 10^-3  found after %d iterations at %g per cent C.L.  "%(  deltar*1000,it_counter, R  )
        return Lmap, deltar, rmax
    
l2=pl.arange(2, 501)
Lcorr1,deltar1,r1= get_delta_r(clc1[2]*1e12 ,lmax=300 , wp=wp, n_steps=2048 )
Lcorr2,deltar2,r2= get_delta_r(clc2[2]*1e12 ,lmax=300 , wp=wp, n_steps=2048 )
In [83]:
fig=pl.figure(figsize=(10 , 5) )

pl.subplot(121)

pl.title('BB')
pl.loglog(l, cl[ 2:lmax+1,3]  ,'k')
pl.xlabel(r'$\ell$', fontsize=16)
pl.ylabel(r'$ C_{\ell}    [ \mu K^2 ]$', fontsize=16)

pl.loglog(l, cl1 [2][2:]  *1e12, alpha=.5, label =r'$f_k =30$ mHz' )
pl.loglog(l, clc1 [2][2:]   *1e12,  label ='w/ mitig. '+r'$f_k =30$ mHz' )

pl.loglog(l, cl2 [2][2:]  *1e12,alpha=.5, label =r'$f_k =150$ mHz')
pl.loglog(l, clc2 [2][2:]   *1e12, label ='w/ mitig. '+r'$f_k =150$ mHz')

pl.loglog(l, pl.ones_like(l)*Nl0,'r', alpha=.5, linewidth=3)
pl.ylim([1e-10,1e-2])
foo=pl.xlim([2,700])
pl.legend() 

pl.subplot(122)

pl.plot([1], [1e-14])
pl.plot(np.linspace(0,.01, 2048), Lcorr1, label=r'$  r=%.1e\pm %.1e   $'%(r1, deltar1 ) )
pl.plot([1], [1e-14])
pl.plot(np.linspace(0,.01, 2048), Lcorr2,  label=r'$  r=%.1e\pm %.1e   $'%(r2, deltar2 ) )

pl.axvline( x=0.57e-3, color='r', alpha=.5, linewidth=3,label='LiteBIRD sys. requirement')
pl.xlim([0., .1e-2])
pl.ylim([0.,1.1 ])
pl.legend(loc='best')
pl.xlabel(r'$r$', fontsize=15)
pl.ylabel(r'$\mathcal{L}(r)$', fontsize=15)
pl.tight_layout()

Fig3 :(left panel) Solid blue and green BB spectra with gain drifts with respectively $30$ and $150$ mHz fknee. Solid orange and red BB spectra post-mitigation for the same cases. (rigth ) Likelihood posterior on $r$ to assess the post-mitigation bias . The thick red line indicate the LiteBird required level from systematics. At the power spectrum level we can observe a similar behaviour as in the maps. The amplitude of BB spectra post mitigation is smaller then the one with the injected gain drifts.

Inputs: Monopole , Dipole, CMB T

We now include the CMB monopole at $2.7 $ K in order to see how the effects of a larger thermal load.

We keep $\sigma_g=0.1$% $f_k=30$ mHz

In [118]:
umap3= hp.read_map(field=range(3) , 
            filename='/Users/peppe/work/satellite_sims/drifting_gain/monopole/uncorr/madam_map_fk0p03.fits', 
                       verbose=False)

fig=pl.figure(figsize=(15,8))
fig.suptitle('gain drifts injected ')

hp.mollview(umap3 [1], unit='K ', title=r'$  Q$', sub=221 ,max=2e-5, min=-2e-5, fig=fig,cmap=cmap  ) 

hp.mollview(umap3[2], unit='K ', title=r'$U$', sub=222 ,max=2e-5, min=-2e-5, fig=fig ,cmap=cmap   )

cmap3= hp.read_map(field=range(3) , 
            filename='/Users/peppe/work/satellite_sims/drifting_gain/monopole/corr/madam_map_fk0p03.fits', 
                       verbose=False)


fig=pl.figure(figsize=(15,8))
fig.suptitle('w/ mitigation')
hp.mollview(cmap3[1], unit='K ', title=r'$  Q$', sub=221 ,max=2e-5, min=-2e-5, fig=fig,cmap=cmap  ) 

hp.mollview(cmap3[2], unit='K ', title=r'$U$', sub=222 ,max=2e-5, min=-2e-5, fig=fig ,cmap=cmap   )

Fig4 (top ) Q and U maps with gain drifts. The presence of monopole causes an increase in the amplitude of fluctuations, (3 orders of magnitude). (bottom) Post-mitigation maps . In this case we don't observe an effective gain with mitigation since the leakage is too noisy for the mitigation to work.

Comment from Reijo:the gain drifts simulated so far are too noisy and not very realistic, the gain is assumed to drift, and it doesn't jump abruplty from a value to another (as noise)

To address this comment I've decreased the amplitude of $\sigma_g=10^{-5} $ and keep the very same $f_{knee} =30$ mHz.

In [121]:
umap3n= hp.read_map(field=range(3) , 
            filename='/Users/peppe/work/satellite_sims/drifting_gain/monopole/uncorr/madam_map_fk0p03_nown.fits', 
                       verbose=False)

fig=pl.figure(figsize=(15,8))
fig.suptitle('gain drifts injected ')

hp.mollview(umap3n [1], unit='K ', title=r'$  Q$', sub=221 ,max=1e-8, min=-1e-8, fig=fig,cmap=cmap  ) 

hp.mollview(umap3n[2], unit='K ', title=r'$U$', sub=222 ,max=1e-8, min=-1e-8, fig=fig ,cmap=cmap   )

cmap3n= hp.read_map(field=range(3) , 
            filename='/Users/peppe/work/satellite_sims/drifting_gain/monopole/corr/madam_map_fk0p03_nown.fits', 
                       verbose=False)


fig=pl.figure(figsize=(15,8))
fig.suptitle('w/ mitigation')
hp.mollview(cmap3n[1], unit='K ', title=r'$  Q$', sub=221 ,max=1e-8, min=-1e-8, fig=fig,cmap=cmap  ) 

hp.mollview(cmap3n[2], unit='K ', title=r'$U$', sub=222 ,max=1e-8, min=-1e-8, fig=fig ,cmap=cmap   )

Fig5 (top ) Q and U maps with gain drifts. As compared with Fig4 here the amplitude of gain fluctuations are 2 orders of magnitude smaller, this is the reason why we observe smaller fluctuations. (bottom) Post-mitigation maps suggest an lower leakage at larger scales , even if they look as noisy as the ones in the top panel.

Let's check how the residual features changes if we increase $f_k= 150$ mHz

In [124]:
umap4= hp.read_map(field=range(3) , 
            filename='/Users/peppe/work/satellite_sims/drifting_gain/monopole/uncorr/madam_map_fk0p15.fits', 
                       verbose=False)

#fig=pl.figure(figsize=(15,8))
#fig.suptitle('gain drifts injected ')

#hp.mollview(umap4 [1], unit='K ', title=r'$  Q$', sub=221 ,max=2e-5, min=-2e-5, fig=fig,cmap=cmap  ) 

#hp.mollview(umap4[2], unit='K ', title=r'$U$', sub=222 ,max=2e-5, min=-2e-5, fig=fig ,cmap=cmap   )

cmap4= hp.read_map(field=range(3) , 
            filename='/Users/peppe/work/satellite_sims/drifting_gain/monopole/corr/madam_map_fk0p15.fits', 
                       verbose=False)
#fig=pl.figure(figsize=(15,8))
#fig.suptitle('w/ mitigation')
#hp.mollview(cmap4[1], unit='K ', title=r'$  Q$', sub=221 ,max=2e-5, min=-2e-5, fig=fig,cmap=cmap  ) 

#hp.mollview(cmap4[2], unit='K ', title=r'$U$', sub=222 ,max=2e-5, min=-2e-5, fig=fig ,cmap=cmap   )
In [126]:
umap4n= hp.read_map(field=range(3) , 
            filename='/Users/peppe/work/satellite_sims/drifting_gain/monopole/uncorr/madam_map_fk0p15_nown.fits', 
                       verbose=False)

fig=pl.figure(figsize=(15,8))
fig.suptitle('gain drifts injected ')

hp.mollview(umap4n [1], unit='K ', title=r'$  Q$', sub=221 ,max=2e-8, min=-2e-8, fig=fig,cmap=cmap  ) 

hp.mollview(umap4n[2], unit='K ', title=r'$U$', sub=222 ,max=2e-8, min=-2e-8, fig=fig ,cmap=cmap   )

cmap4n= hp.read_map(field=range(3) , 
            filename='/Users/peppe/work/satellite_sims/drifting_gain/monopole/corr/madam_map_fk0p15_nown.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
fig.suptitle('w/ mitigation')
hp.mollview(cmap4n[1], unit='K ', title=r'$  Q$', sub=221 ,max=2e-8, min=-2e-8, fig=fig,cmap=cmap  ) 
hp.mollview(cmap4n[2], unit='K ', title=r'$U$', sub=222 ,max=2e-8, min=-2e-8, fig=fig ,cmap=cmap   )

Fig6 (top ) Q and U maps with gain drifts qith a larger fknee (150 mHz). As compared with Fig5, we notice large scale leakage which are reduced in the post-mitigation residuals (bottom) .

In [130]:
cl3=hp.anafast(umap3 , lmax=lmax)
clc3=hp.anafast(cmap3 , lmax=lmax)
cl4=hp.anafast(umap4 , lmax=lmax)
clc4=hp.anafast(cmap4 , lmax=lmax)
In [127]:
cl3n=hp.anafast(umap3n , lmax=lmax)
clc3n=hp.anafast(cmap3n , lmax=lmax)
cl4n=hp.anafast(umap4n , lmax=lmax)
clc4n=hp.anafast(cmap4n , lmax=lmax)
In [136]:
fig=pl.figure(figsize=(10 , 5) )

pl.subplot(121)
pl.tight_layout()
pl.title('BB  '+r'$\sigma_g=10^{-3}$')
pl.loglog(l, cl[ 2:lmax+1,3]  ,'k')
pl.xlabel(r'$\ell$', fontsize=16)
pl.ylabel(r'$ C_{\ell}    [ \mu K^2 ]$', fontsize=16)

pl.loglog(l, cl3 [2][2:]  *1e12, alpha=.5, label =r'$f_k =30$ mHz' )
pl.loglog(l, clc3 [2][2:]   *1e12,  label ='w/ mitig. '+r'$f_k =30$ mHz' )

pl.loglog(l, cl4[2][2:]  *1e12,alpha=.5, label =r'$f_k =150$ mHz')
pl.loglog(l, clc4 [2][2:]   *1e12, label ='w/ mitig. '+r'$f_k =150$ mHz')

pl.loglog(l, pl.ones_like(l)*Nl0,'r', alpha=.5, linewidth=3)
pl.ylim([1e-6,1e1])
foo=pl.xlim([2,500])
pl.legend() 

pl.subplot(122)

pl.title('BB  '+r'$\sigma_g=10^{-5}$')
pl.loglog(l, cl[ 2:lmax+1,3]  ,'k')
pl.xlabel(r'$\ell$', fontsize=16)
pl.ylabel(r'$ C_{\ell}    [ \mu K^2 ]$', fontsize=16)

pl.loglog(l, cl3n [2][2:]  *1e12, alpha=.5, label =r'$f_k =30$ mHz' )
pl.loglog(l, clc3n [2][2:]   *1e12,  label ='w/ mitig. '+r'$f_k =30$ mHz' )

pl.loglog(l, cl4n[2][2:]  *1e12,alpha=.5, label =r'$f_k =150$ mHz')
pl.loglog(l, clc4n [2][2:]   *1e12, label ='w/ mitig. '+r'$f_k =150$ mHz')

pl.loglog(l, pl.ones_like(l)*Nl0,'r', alpha=.5, linewidth=3)
pl.ylim([1e-10,1e-5])
foo=pl.xlim([2,500])
pl.tight_layout()

Fig7 (left) Power spectra of gain drifts T-P leakage from Monopole and Dipole, for two choices of fk 30 (solid blue ) and 150 (solid green) mHz and $\sigma_g=10^{-3}$. (left) As in the left panel, drifts have a lower $\sigma_g=10^{-5}$ . Notice that the mitigation (red and orange solid lines) reduces the leakage at $\ell<10$.